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ABSTRACT 

We present two types of numerical prescriptions that accelerate the radiative transfer 
calculation around point sources within a three-dimensional Cartesian grid by using 
the oct-tree structure for the distribution of radiation sources. In one prescription, 
distant radiation sources are grouped as a bright extended source when the group's 
angular size, ^s, is smaller than a critical value, ^crit, and radiative transfer is solved 
on supermeshes whose angular size is similar to that of the group of sources. The 
supermesh structure is constructed by coarse-graining the mesh structure. With this 
method, the computational time scales with Nm \og{Nm) \og{Ns) where Nm and A^s are 
the number of meshes and that of radiation sources, respectively. While this method 
is very efficient, it inevitably overestimates the optical depth when a group of sources 
acts as an extended powerful radiation source and affects distant meshes. In the other 
prescription, a distant group of sources is treated as a bright point source ignoring the 
spatial extent of the group and the radiative transfer is solved on the meshes rather 
than the supermeshes. This prescription is simply a grid-based version of start by 
Hasegawa & Umemura and yields better results in general with slightly more com- 
putational cost (oc A^m^^ log(A^s)) than the supermesh prescription. Our methods can 
easily be implemented to any grid-based hydrodynamic codes and are well-suited to 
adaptive mesh refinement methods. 
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1 INTRODUCTION 

Radiative transfer (RT) of photons has fundamental impor- 
tance for formation of astronomical objects, such as galax- 
ies, stars, and blackholes. Unfortunately, the nature of RT, 
in which we have to solve the time evolution of the six- 
dimensional phase-space information of photons (three spa- 
tial dimensions, two angular dimensions, and one frequency 
dimension; or equivalent ly three spatial and three momen- 
tum dimensions), makes it difficult to solve RT accurately 
and to couple it with hydr odynamics. To d ate, various RT 
schemes has been proposed (|lliev et al.' 2006), some of which 
are coupled with hydrodynamics (p^liev et al. 200^). A wide 
range of approximation have been used to deal with multi- 
dimensional nature of the transfer equation and they have 
their own pros and cons. 

When radiation sources are embedded in media on 
meshes, RT calculations can be categorised into two types; 
one premises that the source functions are assigned on 
meshes and the other does that radiation sources are treated 
as point sources independent of meshes. In the former type, 
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the RT equations are integrated along long or short charac- 
teristics between meshes. The latter is advantageous when 
the number of the point sources, A^s, is smaller than that of 
the boundary meshes, ^ , where A^m is the total num- 
ber of the meshes. The latter type of the RT schemes is often 
called 'ray-tracing' that we deal with in this paper. 

The most accurate and straight- forward RT scheme is 
the long characteristics method in which all source meshes 
are connected to all other relevant meshes (lAbel et al.lll999l : 
ISokasian et al. I I2OOII : ISusa 2006). This method is however 
very expensive computationally. The computational costs 
scales with in general and with N^^Ns for the transfer 

from point sources. 

The short char act eristic s method (iKunasz Auerl 



1988'- 



iStone et al 



iNakamoto et al.1 l200lh 



[1991; 

reduces 



iMellema et al.l 
the computational 



19981 : 



cost 

by integrating the equation of RT only along lines that 
connect nearby cells. It scales with A^m'^'^ and with NmNs for 
the transfer from point sources. Its known disadvantage is 
the inability to track collimated radiation fields and hence 
the inability to cast sharp shadows owing to the numerical 
diffusion. 

The methods whose computational cost is similar to 
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that of the short characteristics method with smah loss of 
accuracy compared to the long characteristics met hod have 
also been developed ([Razoumov Cardall' '2005' and 'au- 
thentic RT' by N akamoto et al. in llliev et al. 2006, V Adap- 
tive ray tracing (|Abel Wandeltl l2002l) has been w idely 
used for RT around point sources (|Wise Abe]|l201ll ). 

Mote Carlo transport (|Ciardi et alJ 1200 ih is also 
straight forward. The advantage of this approach is that 
comparatively few approximations to the RT equations need 
to be made. The res ulting radiation f ield however inevitably 
becomes noisy (see llliev et al] l2006h due to its stochastic 
nature unless a huge number of photon packets are trans- 
ported. This method is computationally very expensive in 
the optically thick regime. 

The methods, which consider the moments of the RT 
equations and consist in choosing a closure relation to solve 
them, can lead to substantial simplifications that can dras- 
tically speed up the calculations because its computational 
cost scales with ^ A^m- The most common of these meth- 
ods is the flux-limited diffusion, which solves the evolu- 
tion of the first moment and uses a closure relation valid 
in the diffusion limit, which is an isotropic radiative pres- 
sure tensor. The equation is modified with an ad-hoc func- 
tion (the flux limiter) in order to ensure that the radia- 
tive flux is valid in the free-streaming limit. This method 
is very useful in diffu sive regions and ha ve been used to 
stud y accretion discs tohsuga et al. I I2OO5I ) and star forma- 
tion (|Krumholj|2006h . Another method of closing the system 
is the variable Eddington tensor formalism. It gives better 
results than the flux-limited diffusion but are much more 
complex and costly because it requires the local resolution of 
the transfer equation at each timestep. The methods which 
employ t he optically thin varia ble Eddington tensor approx- 
imation (|Gnedin Abell 2001) have been used to study cos- 
mic reionization tanedin Abell200ll : iRicotti et al1l2002l : 
IPetkova & Springel 2009^). A locally evaluated Eddington 
tensor, called the Mi niodel, has also been used to close the 
system ([Gonzalez et al. l2007l) and has app lied to study cos- 
mic reionization (|Aubert TevssierlbOOsf ). The accuracy of 
the moments methods is p roblem-d ependent an d is h ard to 
judge in general situation. IPetkova Springell (|201ll ) have 
developed a method that employs a direct discretisation of 
the RT equation in Boltzmann form with finite angular res- 
olution on moving meshes. This method is advantageous in 
solving problems in which time-dependent solution of the 
RT equation is important. The timestep however has to be 
very short because photons propagate at the speed of light 
unless a reduced speed of light approximation is employed. 

In many astrophysical problems, for example cosmic 
reionization and galaxy form ation, we have to deal with 
numerous radiation sources. IPawlik Schav3 (|2008l ) in- 
troduced source merging procedure in order to avoid com- 
putationally expensive scaling with the number of sources 
and im plemented it on Smooth e d Par ticle Hydrodynamics 
(SPH). iHasegawa Umemural tOld ) utilised the oct-tree 
algorithm ([Barnes Hutlll986l ) in order to accelerate the 
RT around point sources and they coupled the RT with SPH. 
In their method, distant sources from a target gas particle 
are grouped and regarded as a single point source when the 
angular size of the group of the sources is smaller than a 
critical value. Consequently, the effective number of radia- 



tion sources is largely reduced to log(A^s) when there are A^s 
sources. 

The methods we explore in this paper are parallel to this 
approach except that we implement this grouping algorithm 
to grid-based codes. In one of our methods, we introduce su- 
permeshes; a supermesh consists of 8^ meshes and it is char- 
acterised by the mean density of each chemical species of the 
meshes within the supermesh. Solving the RT on superme- 
shes whose angular size is similar to that of the group of 
the sources in question results in further reduction of com- 
putational time in principle. Another approach we take is 
the point source approximation, in which a group of sources 
sufficiently distant from a target mesh is treated as a point 
source. T he latter can be regarded as a grid-based version 
of START ([Hasegawa Umemural [20 id ). 

Unlike gravitational interactions to which the tree- 
algorithm has been widely applied, RT is affected by the 
medium between a source and a target. It is therefore very 
important to test these tree-based approaches in cases where 
an extended group of sources works as a powerful source 
in inhomogeneous medium and affects (e.g. ionizes) distant 
meshes. In this paper, we extensively investigate such cases 
in order to clarify advantages and disadvantages of the meth- 
ods using tree-based algorithm. 

This paper is organised as follows. In section 2, we de- 
scribe the algorithm in detail. In section 3, we present several 
test problems and compare our methods to each other. We 
summarise and discuss the results in section 4. 



2 RADIATIVE TRANSFER WITH 
TREE-ALGORITHM 

In this section, we describe our ray-tracing algorithm that 
we use to solve the steady RT equation for a given frequency. 



(1) 



where lu, r^, and Su are the specific intensity, the optical 
depth, and the source function, respectively. This equation is 
adequate for problems in which the absorption and emission 
coefficients change on timescales much longer than the light 
crossing time. This will always be the case in the volumes we 
will simulate by using our methods. Eqn. ([T]) has a formal 
solution: 



r 

^0 



"dr' 



(2) 



where I^^o is the specific intensity at = and is 
the optical depth at a position along the ray. Throughout 
this paper we employ so-cal led on-the-spot approximation 
([Osterbrock Ferlandl I2OO6I ) in which recombination pho- 
tons are assumed to be absorbed where they were emitted. 
Using the on-the-spot approximation, the formal solution 
given by equation (|2]) is reduced to 



(3) 



To solve this equation numerically, one needs to calculate 
optical depth between each pair of a source and a target 
mesh. The computational cost is hence proportional to the 
number of sources. In the next subsection, we will describe 



the method to decrease the effective number of radiation 
sources by using the oct-tree structure. 

2.1 Source grouping algorithm 

As in iHasegawa Umem ural (|2010l ) , we construct the oct- 
tree structure for the distribution of radiation sources. A cu- 
bic computational domain is hierarchically subdivided into 8 
cubic cells until each cell contains only one radiation source 
or the size of a cell becomes sufficiently small compared 
to that of the computational domain. We call these sub- 
volumes 'tree nodes'. When the side length of the cubic 
computational domain is L, the width of a level / tree node 
is given by w''^^ = L/2K Each tree node records the centre 
of the luminosity of the radiation sources contained in the 
node, 



and the total luminosity, 



(4) 



(5) 



where Vm and Lm indicate the position vector and the lu- 
minosity of a radiation source, respectively, and subscript m 
runs over all sources within the tree node. 

Once we have constructed the tree structure, we loop 
over all meshes. RT from all the radiation sources to each 
target mesh is performed by a simple recursive calculation as 
done in A/'-body calculation. We start at the root node (level 
tree node), which covers entire computational domain. Let 
w be the width of the node currently being processed and D 
the distance between the closest edges of the tree node and 
the target mesh. If the angular size of the node is smaller 
than a fixed value of accuracy parameter, i.e. 



w 
D 



< 



(6) 



then we perform the RT calculation between the group of 
sources in the current node and the target mesh and move 
on to the next node. Otherwise, we examine the child nodes 
(subnodes) recursively. The effective number of sources is 
thus proportional to log(A^s)- In the following subsections 
we will explain how we perform the RT calculation between 
a group of sources and a target mesh. 
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Figure 1. Schematic illustration of the supermesh structure for 
8x8 two-dimensional meshes. In this case, the maximum level, 
^max, is 3 and the meshes themselves can be used as the high- 
est level supermeshes. A level I supermesh contains 2^^^"^^^"^^ 
meshes. For three-dimensional meshes, a level I supermesh con- 
sists of 2'^(^"^a^~^) meshes. 



lo 



Ny |</out>) 



AL 



Figure 2. Plane-parallel radiation with specific intensity Iq enter- 
ing to a supermesh that consists of Nx X Ny meshes. The (i, j)-th 
mesh has the Hi number density, riij. 



supermesh, (/out) (see Fig. [2|. For simplicity, we here only 
consider the absorption by H i atoms and drop the frequency 
dependence. The side length of each mesh is AL and the H I 
number density of the (i, j)-th mesh in the supermesh is riij . 
The mean intensity of the emerging radiation is given by 



(/out) 



^^exp[-crHiA/;-] 



(7) 



2.2 Supermesh approximation 

We first introduce the supermesh approximation. In Fig. [U 
we show a schematic illustration of the supermesh structure. 
If a three-dimensional computational domain is discretised 
by 2^^"^^^ meshes, a level / supermesh consists of 2^*^^"^^^"^^ 
meshes. We can calculate the mean density of each chemical 
species for every supermesh by using the meshes contained in 
it. The meshes can be used as the highest level supermeshes. 
The supermesh structure is resembling to an adaptive mesh 
refinement (AMR) structure and thus this method is well- 
suited to couple with the hydrodynamics by AMR codes. 

Let us consider the case in which plane-parallel radia- 
tion with the specific intensity Iq enters a supermesh that 
consists of Nx x Ny meshes. What we want to know is the 
mean intensity of the ray emerging from the other side of the 



where am is the H i cross-section and A/^ is the H i column 
density of the j-th line, i.e. jVj = J^i^"" rnjAL. 

In our supermesh approximation, we use the mean Hi 
number density, (n) = j ^^ij /{^xNy), to estimate the 
mean intensity of the emerging radiation (/out). Doing this 
introduces some error as we will show below. In order to 
understand the accuracy and the nature of the supermesh 
approximation, we compare the mean intensity of the emerg- 
ing radiation by the supermesh approximation to that calcu- 
lated by using the meshes. We first consider the Taylor series 
expansion of the mean intensity of the emerging radiation 
when we solve the RT on the supermesh: 



(/out) 



7oexp(-(THi(A/'>) 
lo 



(8) 
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where (A/") is the mean H i column density given by 



ALN^{n). 



(9) 



On the other hand, the Taylor series expansion of 
Eqn. © is 



(/out) 



Ny 



E 



-crmN'j + -{crmN'j) 



(10) 



The difference between (/out) and (/out) mean is the second 
order in r. From Eqn. (|7|) and (|8]), the leading error in 

( /out ) mean is 



(/out) ~ (/out) 



out / mean 



/o4^ ((AT^) 



(11) 



Since the variance of the column density, {N"^) — (A/")^, could 
be very large in the inhomogeneous medium, we substan- 
tially overestimate the optical depth if we use Eqn. (|8|). 

We can therefore in principle improve the approxima- 
tion by estimating the variance of the column density. Ac- 
cording to the central limit theorem, the variance of the 
column density for large Nx can be expressed by using the 
variance of the density, if the density, riij, is a sequence of 
independent and identically distributed random variables: 

{Af'} - {Aff = [{n'} - {nf] . (12) 

Using this relation, the mean intensity of the emerging ra- 
diation can be approximated as: 



(/out) 



variance 



/o 



exp {-am{Af)) 



(13) 



The effective column density for a ray segment that inter- 
sects the supermesh is hence 



Met 



In [exp(-aHi(n)/i) + ((n^) - (n)^)] 



(14) 



where h is the length of a ray segment. We however do not 
employ this approximation because Eqn. ()12p is only valid 
for large Nx and Nx always becomes small near the target 
mesh. We thus only use the mean density in our supermesh 
approximation which is described by Eqn. (|8]). We will in- 
vestigate the accuracy of this approximation in Section [S] 

Now we have to determine on which supermeshes we 
perform the RT calculation. We chose to use the lowest level 
supermeshes whose angular size, ^, is equal to or smaller 
than the angular size of a group of the sources, ^s, since 
we assume plane-parallel radiation to construct the approx- 
imation. We define the luminosity- weighted rms projected 
radius as the effective projected size of the group of the 
sourcefl, i.e. if the target mesh is located along the z- 
direction from the centre of the luminosity, the projected 

^ This choice may somewhat underestimate the effective pro- 
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Figure 3. A schematic illustration of the RT calculation using 
supermeshes. The radiation sources in the tree node indicated 
by red square are regarded as a single bright extended source. 
The target mesh is coloured by light blue. The RT is solved on 
the supermeshes at the lowest level, whose angular size is equal 
to or smaller than the angular size of the source group, ^s- The 
supermeshes used are indicated by purple colour and their sizes 
are represented by the orange squares. 



size of the group is defined as 

2 _ X^m {(^m 



(15) 



where x and y are, respectively, the x and y components 
of the position vector of the luminosity centre and the sub- 
script m runs over all sources in the tree node in question. 
Practically, we calculate the following tensor for each tree 
node: 



■ rj)(r„ 



(16) 



where the subscripts i and j, respectively, indicate i-th and 
j-ih components of the position vector, i.e. i and j are either 
X, y, or z; and the subscript m has the same meaning as 
in Eqn. (|15p . By using (0,0) and (1,1) components of the 
tensor X' which is the tensor X in the rotated frame so that 
the target mesh is placed along the z-direction from the 
luminosity centre, we can estimate the angular size of the 
group of the source in the tree node as 



D 



2 
D 



^00 + ^11 



(17) 



where D is the distance between the luminosity centre and 
the closest edge of the target mesh. In Fig. O we illustrate 
the procedure of the RT calculation using the supermeshes. 
The computational cost by this method is expected to scale 

with iV^l0g(iVrn)l0g(iVs). 

jected size as for the case of a disc with a constant surface bright- 
ness. We have confirmed that simulation results are not sensitive 
to such a level of difference (a factor of V2). 
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2.3 Point source approximation 

Here we introduce another way of accelerating the RT calcu- 
lation by using the oct-t ree structure of the distributio n of 
radiation sources. As in iHasegawa Umemural (|2010[ ). we 
treat a group of sources in a tree node which satisfies the 
condition described by Eqn. (|6]) as a bright point source. 
Since we ignore the size of the source group, we solve the RT 
not on the supermeshes but on the meshes. Consequently, 

4 

the computational cost scales with A^m log(As)- Although 
this is slightly more expensive computationally than the su- 
permesh approximation in which the cost is proportional to 
A^m log(A^m) log(A^s), this method may faster than the super- 
mesh approximation for small A^m because we do not have 
to calculate Xqq and X[i in the point source approximatior0- 
Since the surface area of a Stromgren sphere is proportional 
to iV^/^ where N is photoionization rate (see Section [3.ip . 
treating a source group as a point source underestimates the 
surface area of ionized regions. We will explore this effect in 
our tests. 



2.4 Non-equilibrium chemistry 

We solve the non-equilibrium chemistry for e, H I, H II, He i. 
Hell, and Hem implicitly. Note that since we employ the 
on-the-spot approximation, we use 'Case B' recombination 
coefficients to calculate recombination rates of Hii, Hell, 
and He ill throughout this paper. 

Using the optical depth obtained by the methods de- 
scribed in Section 12.21 or 12. 3[ the photoionization rates of 
Hi, He I, and Heii in each mesh are given by 



(18) 



where Fi^a denotes the radiative contribution from a radia- 
tion source (or a group of radiation sources), a, and i = Hi, 
He I, and Hell. The contribution from a point-like radiation 
source, a, is represented by 



47r/ir, 



1 



dzy 



ai{iy)La{iy) exp 



-^A/;-,acr^(zy) 



(19) 

where Ui is the threshold frequency for the i-th species, ai{iy) 
is the cross-section of the i-th species, and r^, Lo,{iy), and 
N'i^a are respectively the distance between the luminosity 
centre and the target mesh, the intrinsic luminosity of the 
radiation source (or the group of the sources), and the col- 
umn density of i-th species. The sum in the exponent runs 
over all three chemical species. When all sources have the 
same spectral shape, i.e. La{h') — Cafi^^), we generate a 
look-up table for each species as a function of column den- 
sities: 



(z/)/(zy) exp 



(20) 



^ It should be noted that, in START (|Hasegawa Umemural 
the computational time scales with Ap log(As), where Ap 
is the number of the SPH particles, by utilising the optical depths 
for SPH p articles in the order of distance from the radiation 
source fsee lSusal[2003 and IHasegawa Umemurall2010l for more 
details). This scaling is better than our point source approxima- 
tion. 



In our case, a look-up table for each chemical species be- 
comes three-dimensional table. We have confirmed that 20 
logarithmic bins for each column density is sufficient. By 
using the look-up tables, the RT calculation is reduced to 
evaluating the column densities. 

Following lAnninos et al ] (I1997I ). we update the densi- 
ties of each chemical species implicitly by using a backward 
difference formula (BDF). The equations to evolve the den- 
sity of each species can be generally written as 



Di (T, rij 



(21) 



where rii — pi / {Aimn) ^ Ai is the atomic mass number of the 
i-th species, and mn is the proton mass. This time i is either 
e. Hi, Hii, Hei, He 11 or Hem. The first term of the right- 
hand side, Ci, is the collective source term responsible for the 
creation of the i-th species. The second term involving Di 
represents the destruction mechanisms for the i-th species 
and are thus proportional to rii. 

Since the timescales for the ionization and recombina- 
tion differ by many orders of magnitude depending on chem- 
ical species, Eqn. ([2T)) is a stiff set of differential equations. In 
numerically solving a stiff set of equations, implicit schemes 
are required unless an unreasonably small timestep is em- 
ployed. As in Anninos et al. (1997.) we adopt a BDF. Dis- 
cretisation of Eqn. ([2T)) yields 



1 + D'+^'At 



(22) 



where all source terms are evaluated at the advanced 
timestep. However, not all source terms can be evaluated 
at the advanced timestep due to the intrinsic nonlinearity 
of Eqn. (|2ip . We hence sequentially update densities of all 
species in the order of increasing ionization states rather 
than updating them simultaneously; We evaluate the source 
terms contributed by the ionization from and recombina- 
tion to the lower states at the advanced timesteps. This 
met hod has been found to be very efficient and a ccurate 
(e.g. lAnninos et"allll997l : lYoshikawa Sasakilliooeh . 

Further improvements in accuracy and stability can be 
made by subcycling the rate solver over a single timestep 
with which the RT is solved. The subcycle timestep, which 
we call the 'chemical timestep', is determined so that the 
maximum fractional change in the electron density is limited 
to 10% per timestep: 



Atchem = 0.1 



(23) 



2.5 Photo-heating and radiative cooling 

Similarly to the photoionization, photo-heating rate for each 
mesh due to the photoionization of the i-th species is given 
by 



(24) 



where 7ii,a indicates the contribution from a radiation 
source (or a group of sources), a, and i = H i. He i, and He 11. 
The total photo- heating rate is defined hy 7i = ^-l-Lini. 
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Table 1. Rates adopted by the code. The hnes are, from top to bottom, reference for: Case B recombination rates (RRB) of Hll, Hell, 
and Hem ; dielectronic recombination rate (DRR) of Hell; coUisional ionization rates (CIR) of Hi, He I, and Hell; Case B recombination 
coohng rates (RCRB) of Hll, Hell, and Helll; dielectronic recombination coohng rate (DRCR) of Hell; coUisional ionization cooling rates 
(CICR) of Hi, He I, and Heii; coUisional excitation cooling rates (CECR) of Hi, Hei, and Heii; bremsstrahlung cooling rate (BCR); 
inverse Compton cooling rate (CCR); photoionization cross-sections (CS) of Hi, He I, and Hell. 



RRB 


DRR 


CIR 


RCRB DRCR 


CICR 


CECR 


BCR 


CCR 


CS 


(4), (5), (4) 


(2) 


(7), (7), (1) 


(4), (5), (4) (3) 


(3), (3), (3) 


(3), (3), (3) 


(4) 


(6) 


(8) 



rPlAbel et al.l (Il997l): ^2) lAldrovandi Peauig^ (119731): ^3) ICe^ (Il992l): (4)lHummed (|l994l ): fS^ lHummer StorevI (|l998l ): (6) 
llkeuchi Ostriked (119861 ): f 7) Ijanev et al.l (ligSTf hfS) lOsterbrock Ferlandl (l2006h 



The contribution from a point-like somce, a, is written as 



47rr2 J V 



(y)Loc(y)(y-yi) exp 



(25) ■ 

As for the photoionization, we generate a look-up table for 
each species when all sources have the identical spectral 
shape. 

We solve the energy equation for each mesh implicitly 



t-\-At t , rt 
U — U -\ 



(26) 



where u and T* = T(n\^u^) are respectively the specific in- 
ternal energy and temperature of the gas and A is the cooling 
function. Although this implicit integration is always stable, 
we need to subcycle the energy solver with Atchem because 
both Ci and Di in Eqn. ([22|) are functions of the tempera- 
ture. We thus perform the rate solver and the energy solver 
alternately. The chemical timestep Atchem is recalculated 
before every subcycle. 



2.6 Chemical reaction and cooling rates 

We try to use the chemical reaction and cooling rates as 
up-to-date as possible. The sources of these rates are sum- 
marised in Table [1] Note that there are notable differences 
in the recombin ation cooling rates between literatures (see 
llliev et aLlliooel ). 



where the second term in the right-hand side prevent the 
timestep from becoming too short when the medium is al- 
most neutral. Our choice for Ce and eni are 0.2 and 0.002, 
respectively. We follow the evolution of the system with the 
minimum of the timestep defined by Eqn. ()27p . i.e. 

At = At.,nain. (28) 

The timestep At is therefore only about twice as long as the 
shortest chemical timestep, Atchem, min- With this timestep, 
we find the solutions typically within 3 to 6 iteration steps. 
While we can of course use a longer timestep, a longer 
timestep requires more iterations and the total number of 
steps becomes similar or even larger than the case we employ 
the timestep defined by Eqn. ([28]). With a longer timestep, 
the solutions sometimes never converge. This timestep is 
in general much shorter than the timestep defined by the 
Courant timestep criterion and therefore we have to subcy- 
cle the hydrodynamical timestep with this timestep when 
we couple the RT with the hydrodynamics. 

When optically thick meshes exist, the solutions con- 
verge very slowly. We thus use smoothed photoionization 
rates, F^, and heating rates, instead of Vi and The 
smoothed rates for the i-ih mesh is calculated by using ad- 
jacent 26 meshes, i.e. 27 meshes in total, with a Gaussian 
kernel of the smoothing length AL. Doing this drastically re- 
duces the number of iterations required to find the solutions. 
This smoothing may introduce the smearing of the I-fronts 
especially when the spacial resolution is poor. While we do 
not find such an effect in our test simulations as we will show 
later, this can be avoided by applyi ng th e smoo thing; only 
to optically thick meshes as done bv lSusal (| 20061 ). 



2.7 Time stepping 

Since the optical depth T(zy) at t+ At depends on densities of 
all species at t + At, we have to solve the static RT equation 
(Eqn. Q), the chemical reactions (Eqn. ([22]) ), and the en- 
ergy equation (Eqn. ([26)) ) iteratively. We iterate these steps 
until the relative difference in the electron number density 
becomes sufficiently small: \ni'^ — rSj^~'^^ \/rSj^^ < e, where 
superscripts indicate the number of iterations and we set e 
to 10~^ throughout this paper. The timestep At, with which 
we solve the RT equation to obtain r.+^* and could 
be much larger than the chemical timestep Atchem, with 
which we subcycle the rate and energy solvers. 

We however choose to employ a timestep that is defined 
by the timescale of the chemical reactions: 



At. 





Tie 


+ em 


nm 


Ce 








Tie 


i 


nm 



(27) 



3 TEST SIMULATIONS 

In this section, we describe the tests we perform. In order 
to evaluate the accuracy of our tree-based RT algorithms, 
problems should involve many sources. Therefore some of 
the tests presented are neither simplest nor cleanest. All test 
problems are solved in three dimensions, with 128^ meshes 
unless otherwise stated. 



3.1 Test 1 — Pure hydrogen isothermal Hii region 
expansion 

The first test is the classical problem of a H II region expan- 
sion in a static, homogeneous, and isothermal gas, which 
consists of only hydrogen, around a single ionizing source. 
This problem has a known analytic solution and is therefore 
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Figure 4. Test 1 - Images of the H l fraction, cut through at the 
mid plane of the simulation box at t = 500 Myr. 



the most widely used test. Note however that since there 
is only a single radiation source, our RT schemes described 
in Sections 12.21 and 12.31 have no difference and both meth- 
ods become the long characteristics method. The aim of this 
test is hence to test our chemical reaction solver and time 
stepping procedure. 

We adopt a monochromatic radiation source that 
steadily emits photons per second, whose frequency is 
the Lyman limit frequency (/iz^l = 13.6 eV). The density 
of the initially neutral gas is nn- Assuming the ionization 
equilibrium, the Stromgren radius is given by 

1/3 



rs = 



3N^ 



(29) 



where is the Case B recombination coefficient. If we as- 
sume that the ionization front (I-front) is infinitely thin, the 
evolution of the I-front radius is analytically given by 

11/3 



n = rs [1 - exp(-t/tr, 



where 



(30) 



(31) 



is the recombination time. 

The analytical solution for the profile of the neu- 
tral and ionized fractions {Xm{r) = nHi(r)/nH and 
Xmi(T) — n\i\\(r) Inn) can also be calculated (e.g. 
lOsterbrock k, Ferlandl fiooil ) from the equation of the ion- 
ization balance at radius r: 

nni{r) 



47rr^ 



^iV^e-"(^^aHi(z^L) 



^ nHii(r) aB(T), 



where 



r(r) = cTHi / niii{r')dr' . 
Jo 

The profile of the neutral fraction is thus given by 



(32) 



(33) 



XHi(r) 



^ \ A O 



H«B J 




r/fs 



Figure 5. Test 1 - The profiles of ionized and neutral fractions. 
The radius is in units of the Stromgren radius. The dot-dot-dot- 
dashed, dotted, dot-dashed, and dashed lines represent simulated 
results at t = 120, 250, 500, and 1000 Myr, respectively. The solid 
line indicates the analytical solution at t = oo given by Eqn. (I34|) . 
The minimum ionized fraction in the numerical results is set by 
the collisional ionization which is not included in the analytical 
solution. 



(34) 



To derive this profile, we ignore the collisional ionization, 
which is included in our simulations. 

The initial physical parameters of this test are the 
same as those of Test 1 in Cosmolo^ic al Radiative Transfer 
Comparison Project (jlliev et al.ll2006l V where the hydrogen 
number density, nn, is 10~^ cm~^, the temperature of the 
isothermal gas is 10^ K, and ionization rate, iV^, is 5 x 10"^^ 
photons s~^. Given these parameters and the recombination 
rate we use, aB (10^ K) = 2.58x10-^2 cm s , the recombi- 
nation time and the Stroemgren radius are tree = 122.6 Myr 
and rs = 5.4 kpc, respectively. 

We employ identical numerical parameters to those in 
llliev et al.l I2OO6I ) : The side length of the simulation box is 
6.6 kpc, initial ionization fraction is set to 1.2 x lO^^, and a 
radiation source is placed at the corner of the box, (0, 0, 0). 
We compare our simulation results to the analytical solution 
given by Eqn. (|34)) which represents the solution at t = 00. 

In Fig.m we show the neutral fraction in the z — 0.5AL 
plane at t = 500 Myr, at which point the I-front is close to 
to the maximum radius, i.e. the Stromgren radius. The H il 
region is nicely spherical, though this is not surprising be- 
cause, with a single source, our method is identical to the 
long characteristics method. In Fig.[5l we show the profiles of 
ionized and neutral fractions at t = 120, 250, 500, and 1000 
Myr. The results asymptotically approach to the analytical 
solution at t = 00. There is a minimum neutral fraction in 
the simulation results, which is set by the collisional ioniza- 
tion that is not included in the analytical solution. 
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Figure 6. Test 2 - Upper panel: Spherically averaged ionized and 
neutral fraction profiles. The dot-dashed, dashed, and solid lines 
indicate indicate the profile at t = 10, 100, and 500 Myr, respec- 
tively. The results from a high-resolution spherically symmetric 
one-dimensional simulation are shown by the dotted lines, which 
almost perfectly overlap with those by the three-dimensional sim- 
ulation. The radius is in units of the Stromgren radius for the uni- 
form isothermal gas with nu = 10""^ cm" -1 andT = 10^ K. Lower 
panel: Spherically averaged temperature profiles. The meaning of 
the lines is the same as in the upper panel. 



3.2 Test 2 — Pure hydrogen Hii region expansion 
with thermal evolution 

Test 2 solves essentially the same problem as Test 1, but 
the ionizing source is assumed to have a 10^ K blackbody 
spectrum and we allow the gas temperature to vary owing 
to heating and cooling processes. The initial gas tempera- 
ture and ionized fraction are set to 10^ K and 1.2 x lO""^, 
respectively. 

In Fig. El we show the neutral and ionized fraction 
profiles (upper panel) and the temperature profiles (lower 
panel) at t = 10, 100, and 500 Myr. We also show the 
results from a high-resolution spherically symmetric one- 
dimensional simulation by the dotted line. For the one- 
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Figure 7. Test 3 - Images of the H l fraction and temperature, cut 
through at the mid plane of the simulation box at t = 30, 100, and 
500 Myr from left to right. The side length of the simulation box 
is 132 kpc in which we randomly distribute 1000 radiation sources 
and 1000 optically thick meshes. Upper two rows show results by 
the supermesh approximation with ^crit =0.5 and lower two rows 



by the point source approximation with 



: 0.5. 



dimensional simulation, we use 1024 meshes for a sphere 
of radius of 1.5 x rs and we do not employ the smoothed 
ionization and heating rates whereas smoothed rates are em- 
ployed in the tree-dimensional simulation. The results by 
the three-dimensional simulation are almost indistinguish- 
able from those obtained by the one-dimensional one. The 
use of the smoothed rates to accelerate the convergence has 
thus no evident side-effects such as smearing of the I- front. 

For this test, our results are most resembling to those 
obtained by rsph for Test 2 in Co smological Radiative 
Transfer Comparison Project (|lliev et al. 2006)0. The agree- 
ment with RSPH is natural because both methods are essen- 
tially the long characteristics method. Small differences are 
probably caused by different adopted rates. 



^ We note that not all codes in Cosmological Radiative Transfer 
Comparison Project were capable of dealing with multifrequency 
RT. 
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Figure 8. Test 3 - Dependence on the accuracy parameter ^crit- 
Upper panels: The volume fractions of the neutral fraction at 
t = 500 Myr. The results by the supermesh approximation are 
presented in the left panel. The solid (black), dotted (red), and 
dashed (blue) lines indicate the results with ^crit = l-O, 0-5, and 
0.0, respectively. The relative difference to the long characteristics 
method (^crit = ^)-> ^5 is also shown. The right panel shows the 
results obtained by the point source approximation. Lower panels: 
The volume fractions of the gas temperature at t = 500 Myr. The 
meaning of the lines are the same as in the upper panels. 



3.3 Test 3 — Multiple radiation sources in a 
clumpy medium 

In order to test the validity of the RT solver based-on the 
source grouping, we have to solve problems that involve 
multiple sources. Moreover, the error in the supermesh ap- 
proximation becomes large when the inhomogeneity of the 
medium is large (see Eqn. (pT]) V In this test, we therefore 
solve the RT from multiple sources in the clumpy medium. 
The side length of the simulation box is 132 kpc. We ran- 
domly select 1000 optically thick meshes whose hydrogen 
number density is nn = 0.2 cm~^ and optical depth at 
the Lyman limit frequency is 4 x 10^ for the mesh size. 
The hydrogen number density of other meshes is set to 
riH = 10 ~^ cm~^. We also randomly distribute 1000 radia- 
tion sources in the simulation box. Each source has a 10^ K 
blackbody spectrum and steadily emits = ^ x 10^^ ion- 
izing photons per second. The initial gas temperature and 
ionization fraction are set to 10^ K and 1.2 x 10""^, respec- 
tively. 

In Fig. [71 we show the neutral fraction and tempera- 
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Figure 9. Test 3 - Relative difference in the temperature, cut 
through at the mid plane of the simulation box at t = 500 Myr. 
This figure compares temperature obtained by the supermesh ap- 
proximation with ^crit = 1 to that by the long characteristics 
method (^crit = 0)- The relative difference in temperature is de- 



fined as At : 



/ rrii supermesh 



^ijlong^^ 2^|long 



ture maps at the mid plane of the simulation box at t = 30, 
100, and 500 Myr. We show the results by the supermesh 
approximation and by the point source approximation with 
^crit =0.5. The results by two methods are virtually iden- 
tical to each other including the shape of shadows by the 
optically thick meshes. 

In order to investigate the dependence on the accuracy 
parameter ^crit, we compare the simulations with ^crit = 1, 
0.5, and 0. In Fig. (8] we show the volume fractions of the 
neutral fraction and the volume fractions of the gas temper- 
ature respectively in the upper panels and lower panels. We 
also show difference in the volume fractions relative to those 
obtained by the long characteristics method (^crit = 0). For 
example, the relative difference in the volume fraction of 
the neutral fraction by the supermesh approximation with 
^crit = X is defined as 



f(^Hi)| 



supermesh 



1/ I long 



f(^Hi)|,, 



(35) 



The volume fractions of the neutral fraction with ^crit = 
1 and 0.5 agree quite well with those by the long character- 
istics method (^crit = 0). The relative differences are typi- 
cally less than 1 % even with ^crit = 1- For a given value 
of the accuracy parameter, the point source approximation 
shows slightly better agreement with the long characteris- 
tics method. On the other hand, agreement in the volume 
fraction of the gas temperature is not as excellent as that for 
the neutral fraction. In particular, both the supermesh and 
point source approximation predict much more low temper- 
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ature gas around 10"^ K. This is because treating a source 
group as a point source underestimates the surface are of 
the ionized regions as we stated in Section [Z3l and the low 
temperature gas is primarily heated by high energy photons 
that permeate beyond the surfaces of highly ionized regions. 
Except for this disagreement for the low temperature gas 
(< 2 X 10^ K), typical difference is less than 10 %. 

To study how serious the deviation from the long char- 
acteristics method at low temperature, we compare the tem- 
perature map obtained by the supermesh approximation 
(^crit = 1), which shows the worst agreement with the long 
characteristics method, and that by the long characteristics 
method in Fig. [9l We find that the temperature difference is 
largest for the low temperature gas with T ^ 10^ K (see also 
Fig.[7j). The difference in temperature is however very small, 
only 10 % at most. We therefore conclude that the results 
with ^crit = 1 are almost converged to the result obtained 
by the long characteristics method. 

This test proves that both tree-based methods produce 
equally good results even with a large value of the accuracy 
parameter, ^crit = 1, in the situation where a local H ii region 
is driven primarily by one or a few sources. This situation is 
resembling to the early stage of cosmic reionization. Only at 
very late stage of the reionization, the H ii regions overlap 
each other and multiple sources become visible each other; 
at this stage, the reionization has largely completed already. 
We thus expect that our tree-based methods, in particular 
the supermesh approximation, are well suited to this type 
of problems. 

3.4 Test 4 — Clustered radiation sources in a 
clumpy medium 

Unlike Test 3, here we explore the problem in which groups 
of sources act like bright extended sources and they ionize 
distant meshes. This would be one of the toughest problems 
for the methods accelerated by source grouping. The side 
length of the simulation box is the same as in Test 3, i.e. 
^box = 132 kpc. In order to construct clustered distribution 
of radiation sources, we put a sphere of radius r — Lbox/4, 
whose centre is randomly placed in the simulation box. We 
uniformly distribute 1000 radiation sources in the sphere. 
We then put a new sphere whose radius is 20% smaller than 
the previous one and again we distribute 1000 sources in the 
sphere. We continue this procedure until we put 10 spheres, 
each of which contains 1000 sources. Consequently, there 
are lO'^ radiation sources in the simulation box. Each source 
has a 10^ K blackbody spectrum and emits Nj = 5 x lO'^^ 
ionizing photons per second. We also randomly select 10^ 
optically thick meshes whose hydrogen number density is 
riH = 0.2 cm""^. The hydrogen number density of the re- 
maining meshes is set to nu — 10~^ cm""^. The initial gas 
temperature and ionization fraction are set to 10^ K and 
1.2 X 10~^, respectively. 

In Fig. IIOI we show the neutral fraction maps, cut 
through at the mid plane of the simulation box. The size and 
shape of the ionized regions by the supermesh approxima- 
tion strongly depend on the value of the accuracy parameter; 
The larger the value is, the smaller the size of the ionized 
regions is. This is due to the very nature of the supermesh 
approximation, which significantly overestimates the optical 
depth when a size of supermesh is large and the variance of 



the Hi density is large (see Eqn. (pTj) and (O). On the 
other hand, the results by the point source approximation 
are relatively insensitive to the value of the accuracy param- 
eter. The size of the ionized regions is almost same between 
^crit = 1 and while small difference is seen in the shapes. 

In Fig. nil we show the volume fractions of the neu- 
tral fraction and gas temperature at t = 100 Myr varying 
the value of the accuracy parameter, ^crit, from 1 to 0. We 
also show the relative difference to the long characteristics 
method (^crit = 0). The volume fraction of the neutral frac- 
tion confirms the dependence of the supermesh approxima- 
tion on the value of the accuracy parameter, i.e. the larger 
the value of ^crit is, the smaller the ionized fraction is. This 
dependence is more evident in the volume fraction of the 
gas temperature. There is more low temperature gas in the 
simulation with a larger value of the accuracy parameter. 
Importantly, the results by the supermesh approximation 
with ^crit = 0.2 still significantly deviate from those by the 
long characteristics methods, and therefore we cannot trust 
the result even with ^crit = 0.2. 

On the other hand, the result by the point source ap- 
proximation with ^crit = 1 shows an excellent agreement 
with that with the long characteristics method, in spite of 
the fact that this approximation ignores the spatial extent 
of source groups. This result proves that the point source 
approximation is very efficient and accurate for this type of 
problems. 

The relative difference to the long characteristics 
method indicates that both approximations overestimates 
the volume fraction of the almost fully-ionized gas (^fj j — 
2 X 10~^). This ionized fraction corresponds to the central 
regions of each source spheres. The volume of these regions 
are however very small and the neutral fraction is very low 
anyway; this overestimation of the ionization fraction at the 
central regions of the source spheres does not affect the evo- 
lution of the whole simulation box. In fact, by the point 
source approximation, the relative difference to the long 
characteristics method in the volume fraction of the neu- 
tral fraction is typically 1 % and ~ 10 % at most except for 
the highly ionized gas with -^fjj < 10~^. 

Even by the point source approximation, the relative 
difference in the volume fraction of the gas temperature to 
the long characteristics method is rather large for the low 
temperature gas. The gas temperature however agrees very 
well with that by the long characteristics method just as we 
showed for Test 3. Except for the low temperature gas, the 
typical difference is ~ 10 %. Interestingly, decreasing the 
value of the accuracy parameter in the point source approx- 
imation from 1 to 0.2 does not improve the agreement with 
the long characteristics method very much in spite of the 
fact that the simulations with a smaller value of the accu- 
racy parameter is much more computationally expensive as 
we will show in the next subsection. Since the point source 
approximation with ^crit = 1 seems to be sufficiently ac- 
curate, we expect that this approximation with ^crit = 0.5 
would be a safe choice for most types of problems. 

3.5 Code performance 

We here investigate how the computation time scales with 
the number of meshes and that of the sources. For this pur- 
pose, we measure the wall-clock time taken for one step of 
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Figure 10. Test 4 - Hi fraction maps, cut through the simulation box at coordinate z = 63. SAL = 65.5 kpc at time t = 100 Myr. From 
left to right, the values of the accuracy parameters are ^crit = l-O, 0-5, 0.2, and 0.0 respectively. The upper panels show results by the 
supermesh approximation and the lower panels display those by the point source approximation. 



the RT calculation. The computation time for solving chem- 
istry etc. is not included. We use 8 cores of 2.13 GHz Xeon 
E5506 processors for these simulations. 

In order to study the scaling with the number of the 
meshes, we randomly place 1000 radiation sources in the 
simulation box. Each source and the simulation box is the 
same as used in Test 1 except that there are 1000 sources and 
we vary the number of the meshes. We show the result in the 
upper panel of Fig. [121 We find that the supermesh approxi- 
mation is slightly faster than the point source approximation 
for a given set of A^m and ^crit • The computation time by the 
point source approximation is a slightly steeper function of 
the number of the meshes than that by the supermesh ap- 
proximation. The computation time by the point source ap- 
proximation scales with A^m'^^ as expected. The scaling of the 
computation time by the supermesh approximation is some- 
where between oc A^m log(A^m) and oc N^l^ . Since the RT is 
solved on the supermeshes whose angular size is similar to 
the angular size of the source group, ^s, which can be much 
smaller than ^crit, the computation time becomes steeper 
function of A^m than the expected scaling, oc A^m log(A^m)- 

In the lower panel of Fig. I12[ we plot the computation 
time as a function of the number of the sources. The number 
of the meshes is fixed to 128"^. The computation time scales 
with log(A^s) for A^s > 1000 in all cases. This result proves 
that the tree-based source grouping is quite efficient to deal 
with a large number of radiation sources. For a given set of 
A^s and ^crit, a simulation by the supermesh approximation 
is always faster than that by the point source approxima- 
tion. It should be however noted that even with the same 
value of the accuracy parameter, simulations by the point 



source approximation are sometimes much more accurate 
than those by the supermesh approximation as we showed 
by Test 4. 



4 SUMMARY AND DISCUSSION 

We have presented a code to solve radiative transfer around 
point sources within a three-dimensional Cartesian grid, ar- 
got, which accelerates the RT calculation by utilising the 
oct-tree structure in order to reduce the effective number 
of radiation sources. We have explored two methods: one 
is the supermesh approximation and the other is the point 
source approximation. In both methods, sources in a tree 
node whose angular size is smaller than the accuracy pa- 
rameter ^crit are treated as a single bright source. As a re- 
sult, computation time only scales with log(A^s). The main 
difference between these two method is that while the for- 
mer takes the spatial extent of a source group into account, 
the latter ignores the size of the source group and treat it 
as a point source. In the supermesh approximation, the RT 
is solved using supermeshes whose angular size is similar to 
the angular size of the source group in question. Doing this 
results in the further acceleration of the RT calculation. 

One might thus see that the supermesh approximation 
is superior to the point source approximation. We have how- 
ever shown that the point source approximation is always 
equally or more accurate than the supermesh approxima- 
tion for a given value of the accuracy parameter. This is 
because RT in a inhomogeneous medium on a supermesh 
inevitably overestimates the optical depth. This approxima- 
tion can be in principle improved by including higher order 
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Figure 11. Test 4 - Dependence on the accuracy parameter ^crit- 
Results at t = 100 Myr are presented. The results are displayed 
in the same manner as Fig. [8] The volume fractions for the sim- 
ulation with ^crit = I5 0-5, 0.2, and are indicated by the solid 
(black), dotted (red), dot-dashed (green), and dashed (blue) lines, 
respectively. 



moments, such as variance, although we do not take such an 
approach. This method hence only applicable to the prob- 
lems in which a local H 11 region is driven primarily by one or 
a few sources such as Test 3 in this paper. When one applies 
the supermesh method to the simulation of cosmic reioniza- 
tion, it could be comb ined with the 4ocal clump ing factor' 
approach proposed bv iRaicevic TheunsI diml ). although 
exploring such a method is beyond the scope of this paper. 

The point source approxi mation, which can be regarde d 
as a mesh version of start ([Hasegawa Umemural]201Ql ). 
produces sufficiently accurate results with ^crit = 1 for all 
test simulations presented in this paper. This approximation 
requires slightly more computational cost than the super- 
mesh approximation and it scales with A^m'^^ log(A^s)- The 
performance can be improved if we choose the angular res- 
olution so that at least one ray from a radiation source (or 
a group of sources) crosses all target meshes instead of solv- 
ing RT to all target meshes. Doing this reduces the total 
number of rays from oc Nm to oc . Such an algorithm 
has b een applied for RT from point sources ([Yaiima et al.l 
and can be extended to our tree-based algorithm. The 
expected scaling is A^m log A^s , which is even faster than the 
supermesh approximation and the same scaling by start. 

For parallel implementation, if the entire meshes and 
sources can fit into memory of one computer node, paral- 
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Figure 12. Computation time taken for one step of the RT cal- 
culation. Upper panel'. Computation time as a function of the 
number of meshes, A^m- The number of radiation sources, A^s is 
fixed to 1000. The solid (red) and dashed (green) lines show the 
results by the supermesh approximation with ^crit = 1 and 0.5, 
respectively. The dot-dashed (blue) and dotted (light blue) lines 
indicate the point source approximation with ^crit = 1 and 0.5, 
respectively. The thin dot-dot-dot-dashed lines show the scaling 
with A^m'^^ and A^mlog(A^m)- Lower panel: Same as the upper 
panel but the number of the source, A^s, is varied. The number 
of meshes, A^m is fixed to 128^. The thin dot-dot-dot-dashed line 
indicates the scaling with log(A^s)- 



lelisation via angle decomposition is preferable to volume 
decomposition. We implement the angle decomposition by 
using both mpi and OpenMP. If a simulation size becomes too 
large to fit into the memory of one computer node, we have 
to employ the volume decomposition. The volume decom- 
ositi on for RT around point sources was introduce by[^ 



i 



20061 ) and the algorithm can be applied to our methods. 
We leave the volume decomposition to future work. 

The method presented in this paper can be easily com- 
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bined with any ^rid-ba sed hydrodynamic c ode, even with 
codes based on A MR (iFrvxell et all l2QQ0l : iTevssierl l20Q2l : 
lO'Shea et aDl2004l ) and will be useful for various astrophysi- 
cal problems in which a large number of radiation sources are 
required such as cosmic reionization and galaxy formation. 
We will apply our code for these issues in a forth coming 
paper. 
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